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The diffusion and flow of amorphous materials, such as glasses and granular materials, has resisted a simple 
microscopic description, analogous to defect theories for crystals. Early models were based on either gas-like 
inelastic collisions or crystal-like vacancy diffusion, but here we propose a cooperative mechanism for dense 
random-packing dynamics, based on diffusing "spots" of interstitial free volume. Simulations with the Spot 
Model can efficiently generate realistic flowing packings, and yet the model is simple enough for mathematical 
analysis. Starting from a non-local stochastic differential equation, we derive continuum equations for tracer 
diffusion, given the dynamics of free volume (spots). Throughout the paper, we apply the model to granular 
drainage in a silo, and we also briefly discuss glassy relaxation. We conclude by discussing the prospects of 
spot-based multiscale modeling and simulation of amorphous materials. 



Dedicated to Prager Medalist, Salvatore Torquato. 



1. Introduction 

Professor Torquato has made pioneering con- 
tributions to the characterization of random 
packings and their rel ation to p r opert ies of het- 
erogeneous materials l)Torauatol [2.003) . His re- 
cent work rejects the classical notion of "random 
close packing" of hard spheres and replaces it 
with the more precise co ncept of a "max i mally 
random jammed state" llTorauato et all [2 000; 



ITorauato and Stillingerl l200l[ iKaris^deT^d 



120021 iDonev et alll2004bl). In the se studies and 
others l|0'Hern et all l200l 12003)) . however, the 
focus is on the statistical geometry of static pack- 
ings, and not on the dynamics of nearly jammed 
packings in flowing amorphous materials. 

Dense random-packing dynamics is at the heart 
of condensed matter physics, and yet it remains 
not fully understood at the microscopic level. 
This is in contrast to dilute random systems 
(gases), where Boltzmann's kinetic theory pro- 
vides a successful statistical description, based on 
the hypothesis of randomizing collisions for indi- 
vidual particles. The same single-particle theory 
can also be applied to molecular liquids at typ- 
ical experimental time and length scales, where 
kinetic energy is able to fully disrupt local pack- 
ings. The difficulty arises in describing liquids 
at very small (atomic) length and time scales, 



over which the trajectories of neighboring par- 
ticles a re strongly correlated - the so -called cage 
effect (jHansen and McDonald 11986). 

This difficulty is extended to much larger 
length and time scale s in glassy relax- 
ation lAn^el^t_alJJyOO0|) and dense granular 
flow l|Jaeger et all ll9961. where kinetic energy is 
insufficient to easily tear a particle away from its 
cage of neighbors. As a result, one must somehow 
describe the cooperative motion of all particles 
at once. In dense ordered materials (crystals), 
cooperative relaxation and plastic flow are medi- 
ated by defects, such as interstitials, vacancies, 
and dislocations, but it is not clear how to define 
"defects" for homogeneous disordered materials. 

The challenge in describing random-packing 
dynamics is related to the concept of "hyper- 
unif orm" point distributions , recen tly introduced 
by ITorauato and Stillingerl l|2003|) . In a dilute 
gas, particles undergo independent random walks 
and thus have the "uniform" distribution of a 
Poisson process ijHadiiconstantinou et alll2003|) . 
where the variance of the number of particles, N, 
scales with the volume, V: Var(iV) = (N) = pV, 
where p is the mean density. In a condensed 
phase, the particle distribution must be "hyper- 
uniform" , with much smaller fluctuations, pro- 
portional to the surface area: Var(iV) cx V^ 1 ^ d , 
where d = 3 is the spatial dimension, so it is clear 
that particles cannot fluctuate independently. In 
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Figure 1. Eyring's mechanism for flow in viscous 
liquids: A single particle jumps from one available 
"cage" to another by exchanging with a "void" 
moving in the opposite direction. 



a crystal, hyper-uniformity is a property of the 
ideal lattice, which is preserved during diffusion 
and plastic flow by the motion of isolated defects. 
For dense disordered materials, however, no sim- 
ple flow mechanism has been identified, which 
p reserves hyper- uniformity. 

lEvrind l|l936f) was perhaps the first to sug- 
gest a microscopic mechanism for viscous flow in 
liquids, analogous to vacancy diffusion in crys- 
tals. He proposed that the packing evolves when 
individual particles jump into available cages, 
thus displacing pre-existing voids, as shown in 
Fig. Much later, the same hypothesis was put 
forth i ndependently in theor i es of t he glass tran- 
sition l(Cohen and Tiirnbull 119591) . shear flow 
in metallic glasses |SpaenenL Il977l). granular 
drainage from a silo (Litwiniszyn, 1958; Mullins, 
119721) . and compaction in vib rated granular ma- 
terials (|Boutreux et all Il998l) . although it is not 
clear that all of these authors intended for the 
model to be taken literally at the microscopic 
level. Since particles and voids simply switch 
places, it seems the Void Model can only be simu- 
lated on a single, fixed configuration of particles, 
but real flows are clearly not constrained in this 
way. This difficu l ty is apparent in the work of 
ICaram and Honel lll99l|) . who neglected random 
packings and simulated the Void Model on a lat- 



tice in an attempt to describe granular drainage. 

By now, free volume theories of amor- 
phous materials (based on voids) have fallen 
from favo r, and expe ri ments on glassy re- 
laxation llWeeks et all l2000l) and granular 
drainage flChoTeT all 12004^) have demonstrated 
that packing rearrangements are highly cooper- 
ative and not due to single-particle hops. In a 
recent theory of granular chute flows down in- 
clined planes, lErtas and Halsevl l|2002j) have pos- 
tulated the existence of coherent rotations, called 
"granular eddies" to motivate a continuum theory 
of the mean flow. Although the theory success- 
fully predicts Bagnold rheology and the c ritica l 
layer thickness for flow. lLandrv and Grestl |2003) 
have failed to find any evidence for granular ed- 
dies in discrete-elem ent simulations of chute flow 
I Landrv et allEoO? 



In glass theory, 
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troduced the concept of regions of coopera- 
tive relaxation, whose sizes diverge at the glass 
transition. Modern statistical mechanical ap- 
proaches are based on mode-coupling theory 
ijGotze and Siogre"rl ^)92) , which accurately pre- 
dicts density corr elation functions in simple liq- 
uids llKobl Il997ft . although a clear microscopic 
mechanism, which could be used in a particle 
simulation, has not really emerged. Coopera- 
tive rearrangements have also long bee n recog- 
nized in the literature on sheared glasses. lOrowanl 
(1952) was perhaps the first to postulate lo- 
calized shear transforma tions in regio ns of en- 
hanced atomic disorder. lArgonl Jl979l) later de- 
veloped the idea of "intense shear transforma- 
tions" at low temperature, which underlies the 
stochastic model of "localized inelas tic transfor- 
mations" l|Bulatov and Argorl Il994j) . A similar 
notion of "shear trans formation zones" (STZ) has 
also been invoked by iFalk and Langerl l)l998l) in 
a continuum theory of shear response, which has 
recently been extended to account for fr ee- volume 
creatio n and annihilation in glasses jLemaitrel 
2002b) and granular materials l)Lemaitrel 120 02a) . 
This phenomenology seems to capture many uni- 
versal features of amorphous materials, although 
the microscopic picture of "+" and "— " STZ 
states remains vague. 

In this paper, we propose a simple model for 
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the kinematics of dense random packings. In 
section |21 we introduce a general mechanism for 
structural rearrangements based on the concept 
of a diffusing "spot" of free volume. In section|31 
we apply the Spot Model to granular drainage. 
In section we analyze the diffusion of a tracer 
particle via a non-local, nonlinear stochastic dif- 
ferential equation, in the limit of an ideal gas of 
spots. In sectional we derive equations for tracer 
diffusion in granular drainage, which depend on 
the density, drift, and diffusivity of spots (or free 
volume). We close by discussing possible applica- 
tions to glasses in sectionEland spot-based multi- 
scale modeling and simulation of amorphous ma- 
terials in section [3 

2. The Spot Model 

2.1. Motivation 

Our intuition tells us that a particle in a dense 
random packing must move together with its 
nearest neighbors over short distances, followed 
by gradual cage breaking at longer distances. In 
simple liquids, this transition occurs at the molec- 
ular scale (< nm) over very short times (< ps) 
compared to typical experimental scales. In su- 
percooled liquids and glasses, the time scale for 
structural relaxation effectively dive rges and is re- 
placed by slow, power-law decay (lAngell et all 
1200ft iKohl Il997t lHansen and McDonald! Il98(j) . 
although the length scale for cooperat i ve mo tion 
remains relatively small llWeeks et a.1.1 l200Ctl . In 
granular drainage, cage breaking occurs slowly, 
over time scales comparable to the exit time from 
the silo, so that cooperative motion is important 
throughout the sy stem at the macroscopic scale 
llChoi et all 120041) . 

Another curious feature of granular drainage is 
the importance of geometry: all fluctuations in a 
dense flow seem to have a universal dependence 
on th e distance droppe d for a wide range of flow 
rates ijChoi et alll2004|) . In a sense, therefore, in- 
creasing the flow speed in this regime is like fast- 
forwarding a film, passing through the same se- 
quence of configurations, only more quickly. The 
only existing theory to predict this property (as 
well as the mean flow profile in silo drainage) is 
the Void Model, since increasing the flow speed 



simply increases the injection rate of voids, but 
not their geometrical trajectories. However, the 
model incorrectly predicts cage breaking and mix- 
ing at the scale of individual particles. This 
"paradox of granular diffusion" is a fruitful start- 
ing poin t for a new model of random-packing dy- 
namics IjBazant et all 12005') . 

2.2. General formulation 

Let us suppose that the cage effect gives rise 
to spatial correlations in particle velocities, with 
correlation coefficient, C(r), for two particles sep- 
arated by r. More generally, when there is bro- 
ken symmetry, e.g. due to gravity in granu- 
lar drainage, there is a correlation coefficient, 
Cp' 3 (ri,r 2 ), for the a velocity component of a 
particle at ri and the f3 velocity component of a 
particle at r 2 . Perhaps the simplest way to en- 
code this information in a microscopic model is to 
imagine that particles move cooperatively in re- 
sponse to some extended entity - a spot - which 
causes a particle at r p to move by 

AR P = -w(r p ,r s ) AR S (1) 

when it moves by AR S near r s . Although the 
spot is not a "defect" per se, like a dislocation 
in an ordered packing, it is a collective excitation 
which allows a random packing to rearrange. 

In principle, the spot influence, w, could be 
a matrix causing a smooth distribution of local 
translation and rotation about the spot center r s , 
but a reasonable first approximation is that it is 
simply a collective translation the opposite direc- 
tion from the spot displacement, as illustrated in 
Fig. In this case, w is a scalar function, whose 
shape is roughly that of the velocity correlation 
function. More precisely, under some simple as- 
sumptions, we show below that C(r) is the over- 
lap integral of two spot influence functions sepa- 
rated by r. Since the spot influence is related to 
the cage effect, we expect that w and C will decay 
quickly with distance, for separations larger than 
a few particle diameters. Due to the local statisti- 
cal regularity of dense random packings, we might 
expect a spot to retain its "shape" as it moves, in 
which case w depend only on the separation vec- 
tor, y 1 ^ — rg\ although this assumption might 
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Figure 2. The spot mechanism for cooperative 
diffusion: A group of neighboring particles makes 
small correlated displacements in response to a 
diffusing "spot" of excess interstitial volume. 



Figure 3. The trails of spots correspond to tran- 
sient, reptating chains of particles. 



need to be relaxed in regions of large gradients in 
density or velocity. 

Physically, what is a spot? Since particles move 
collectively in one direction, a spot must corre- 
spond to some amount of free interstitial volume 
(or "missing particles" ) moving in the other di- 
rection. If particles are distributed with number 
density p P (v p ) and a spot at r s carries a typical 
volume, V s (r s ), then an approximate statement 
of volume conservation is 

V S AR S = - J dr p p(r p )w(r p , r s )AR p (r p ). (2) 

For particles distributed uniformly at volume 
fraction, </>, this reduces to a simple expression 
for the local volume carried by a spot, 

V s (r s )=(f> J dr p w(r p ,r s ), (3) 

which can thus be interpreted as a measure of the 
spot's "total influence". 

A very simple, spot-based Monte Carlo simula- 
tion proceeds as follows. Given a distribution of 
(passive) particles and (active) spots, the random 
displacement, AR^, of the jth spot centered at 



r s induces a random displacement, ARj/, of the 
ith particle centered at r p : 

ARW = -w(r® , t& + AR[ j) ) AR< J > (4) 

Each spot undergoes an independent random 
walk, with an appropriate drift and diffusivity for 
free volume, which leaves in its trail a thick chain 
of particles reptating in the opposite direction, as 
shown in Fig. 03 In Eq. we choose to center 
the spot influence on the end of its small displace- 
ment, but it is also r easonable to use the midpoint 
of the displacement l|Rvcroft et alll2005j) . In the 
infinitesimal limit (see below), these choices are 
analogous to differ ent de finitions of stochastic dif- 
ferentials l|RiskenLll99^ . 

In principle, the drift velocity, diffusivity, and 
influence function of spots could depend on lo- 
cal variables, such as stress, and temperature 
(or rather, some suitable microscopic quantities 
related to contact forces and velocities, respec- 
tively). Spots could also interact with each other, 
undergo creation and annihilation, and possess a 
statistical distribution of sizes (or influence func- 
tions). The simplest kinematic assumption, how- 
ever, which captures the basic physics of the 
cage effect, is that spots are identical and main- 
tain their properties while undergoing indepen- 
dent (non-interacting) random walks. In particu- 
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lar, the constant influence function, w(\r p — r s |), 
is chosen to be translationally invariant in space 
and time. It turns out that this model allows 
rather realistic multiscale simulations, while re- 
maining analytically tractable. 



2.3. Multiscale simulation 

The simple spot mechanism above gives a rea- 
sonable description of tracer diffusion and slow 
cage breaking in random packings, but it does 
not strictly enforce packing constraints (or, more 
generally, inter-particle forces). As such, parti- 
cles perform independent random walks in the 
long-time limit, which eventually leads to uniform 
density fluctuations with Poisson statistics. For 
a complete microscopic model, we must somehow 
preserve hyper-uniform packings. 

This may be accomplished by adding a relax- 
ation step to the spot-induced displacements, as 
shown in Fig. First in (a) , a spot displacement 
causes a simple correlated displacement, as de- 
scribed above, e.g. using Eq. |@J with some short- 
ranged choice of w(r) with a finite cutoff. Next 
in (b), the affected particles and a shell of their 
nearest neighbors are allowed to relax under ap- 
propriate inter-particle forces, with more distant 
particles held fixed. For simulations of (nearly) 
hard grains, the most important forces come from 
a soft-core repulsion, which pushes particles apart 
only if particle begin to overlap. 

Although it is not obvious a priori, the net 
spot-induced cooperative displacements, shown 
in Fig. Efc), easily produce very realistic flow- 
ing packings, whil e preserving the phy sical pic- 
ture of the model i Rvcroft et allEooiK . In prac- 
tice, the correlated nature and small size of the 
spot-induced block displacements results in very 
small and infrequent particle overlaps, only near 
the edges of the spot, where some shear occurs 
with the background packing. As a result, it 
seems the details of the relaxation are not very 
important, although this issue merits further in- 
vestigation. In any case, the algorithm is inter- 
esting in its own right as a method of multiscale 
modeling, since it combines a macroscopic sim- 
ulation of simple extended objects (spots) with 
localized, microscopic simulations of particles. 



3. Application to Granular Drainage 

3.1. Spot parameters 

The classical Kinematic Model for the mea n 
velocity in granular drainage l|Neddermanlll99l|) . 
which compares fairly well with experiments 
dTiiziin and NeddermanL Il979t ISamadani et all 
ll999ilChoi et all 1200511 . postulates that the mean 
downward velocity, v, satisfies a linear diffusion 
equation, 



— = b V ± v 

oz 



(5) 



where the vertical distance z plays the role of 
"time" and the horizontal dimensions (with gra- 
dient, VjJ play the role of "space". The mi- 
croscopic justification for E q. @ is t he Void 
Model of iLitwiniszvnl <|l95Sl |l963) and iMullinsI 
lll972l I1974IT where particle-sized voids perform 
directed random walks upward from the orifice. 
As discussed above, this microscopic mechanism 
is firml y contradicted by particle-tracking experi- 
ments l|Choi et all 1200^1 . but, as shown below, a 
similar macroscopic flow equation can be derived 
from the Spot Model, where spot diffuse upward 
with a (horizontal) diffusion length, 



Var(Ax s 
2d h Az s 



(0) 



where Ax s is the random horizontal displacement 
of a spot as it rises by Az s and dh = 2 is the 
horizontal dimension. A typical value for 3 mm 
glass beads is b « 1.3d, where d is the particle 
diameter. 

The shape of the spot influence function can 
be inferred from measureme nts of spatial velocity 
correlations in experiments ijBazant et all [2005) 
or simul ations by the discre te-element method 
(DEM) l|Rvcroft et all l2005|) . The simplest as- 
sumption is a uniform spherical influence with a 
finite cutoff, 



w(r) 



w 




r < d s /2 
r > d s /2 



(7) 



where experiments and simulations find d s w 5d. 
This is consistent with our interpretation of the 
spot mechanism in terms of the cage effect, where 
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Figure 4. Multiscale simulation of densely packed (nearly) hard spheres with the Spot Model: (a) A 
block of neighboring particles translates opposite to the displacement of a spot of free volume; (b) the 
block and a shell of neighbors are allowed to relax under soft-core repulsive forces; (c) the net cooperative 
rearrangement combines these two steps. (Particle displacements are greatly exaggerated for clarity.) 



a particle moves with its nearest neighbors. The 
typical number of particles affected by a spot, 



N = 



(8) 



is thus N w 72, for <j> « 0.58. 

For a uniform spot, the condition of volume 
conservation, Eq. ©, reads 



V s (Ax,, Az p ) = -NV P (Ax p , Az p ). 



(9) 



Using Eq. QJ, this provides an expression for the 
spot influence, 



Vs 



Acj) 



(10) 



in terms of A(j>, the change in local volume frac- 
tion due to the presence of a single spot. In 
discrete-element simulations of granu lar spheres 
in silo drainage l|Rvcroft et all lioOsT) . the local 
volume fraction varies in the range, (f> = 0.565 — 
0.605, within the rough bounds of jammin g 



0.63 llTorauato et al 



l2000t iKansal et all 1200 



lO'Hern et all 120021 200.C and random loose 
packing, cf) = 0.55 l|Onoda and Linigeri [1.990). 
If we attribute A(j)/<f> = 1% to a single spot, 
then we find w w 0.017, but, if many spots, say 



N s — 10, can overlap, then this estimate is re- 
duced by l/N s . We thus expect, w = 10~ 3 -10~ 2 , 
which can be tested against experiments. 

We can also use Eqs. (I9*)l- H10[l to infer w from a 
measurement of the horizontal particle diffusion 
length: 



Var(Axp) _ w 2 Var(Ax s 



P 2d h \Az p \ 

Discr ete-element 
120051 and 



2dhwAz s 



w b 



(11) 



simulations l|Rvcroft et all 
particle-tracking experiments 



llChoi et all 12001 for similar flows yield 
w = b p /b = 0.00286d/1.14d = 0.00250 and 
w = d/bPe x = d/(1.3d)(321) = 0.0024, respec- 
tively (where Pe x is a Peclet number). These 
values are consistent with the predictions of the 
model. 

A spot's influence is related to the free volume 
it carries by Eq. © . In the case of a uniform spot 
of diameter, d s , its total free volume is given by 



V s 



ir^wd^ 
6 



(12) 



which is related to the particle volume, V p , by 



Vn 



w<p 



= wN. 



(13) 
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For dense granular drainage, the typical values 
N = 72 and w = 0.0025 imply V s = 0.18V P , so 
a spot carries only around one fifth of a parti- 
cle volume, spread out over a region of roughly 
five particle diameters. This delocalized picture 
of free volume diffusion is radically different from 
the classical Void Model, which is the (unphysi- 
cal) limit where V s — V p and N = w = 1. 

3.2. Spot simulations 

As described above, the Spot Model can be suc- 
cessfully calibrated for granular drainage to en- 
able very simple and efficient simulations. First, 
we consider the basic mechanism in Fig. El and 
Eq. applied to granular drainage in a quasi- 
two-dimensional silo with a narrow opening, far 
from the side walls, a s in the experiments o f 
ISamadani et~aTl l|l999l) and IChoi et all (|2004h . 
For simplicity, we simulate the model in two di- 
mensions (dh = 1) using uniform spots with the 
parameters determined from experiments above 
(w = 0.0024, d s = 5d,b= 1.3d). The simulation 
begins with a random packing of identical disks, 
colored with horizontal stripes (10d thick) to aid 
in visualizing the subsequent evolution. 

A snapshot of the spot simulation at a later 
time is shown in Fig. |5| a). For comparison, a 
simulation of the same situation with the Void 
Model on a two - dimen sional lattice, following 
ICaram and Honel l|l99ll) . is shown in Fig. Efb), 
along with the central slice of a three-dimensional 
DEM simulation (15 d thick) in Fig. ^Jd), which 
is very similar to the experiment. Although the 
mean flow profile is similar in the spot and void 
simulations and reasonably close to experiment, 
the void simulation displays far too much diffu- 
sion and cage breaking, since the initial horizontal 
stripes are completely mixed down to the single 
particle level inside the flow region. In contrast, 
the interfaces between the colored layers remain 
fairly sharp in the spot simulation, as in experi- 
ments and the DEM simulations. Tracer diffusion 
is described fairly well by the spot simulation, al- 
though the lack of packing constraints eventually 
leads to the loss of hyper-uniformity, as particles 
begin to overlap and open gaps in the lower side 
regions of highest shear. For more details, see 
iBazant et al"f[|2005h . 



Next, we consider the multiscale spot algorithm 
in Fig. 01 with an internal relaxation step, applied 
to a three-dimensional drainage simulation, start- 
ing from the same initial condition and geometry 
as the DEM simulation in Fig. E^d). The sim- 
plest possible relaxation scheme is to push any 
pair of overlapping particles in a relaxation zone 
(a sphere of diameter, d s +4d) apart by a displace- 
ment, ct(d — r) proportional to the overlap, d — r, 
while keeping particles fixed outside a sphere of 
diameter, d s + 2d. A spot simulation with similar 
parameters as above, including such a relaxation 
step with a = 0.8, is shown in Fig. Efc). The 
rate of introducing spots at the orifice and their 
upward drift velocity has also been calibrated for 
comparison to the DEM simulation, at the same 
instant in time. 

Clearly, the simple multiscale relaxation step 
is able to preserve realistic random packings, and 
in many ways the spot simulation in Fig. Efc) 
is indistinguishable from the much more com- 
putationally demanding simulation in Fig. |Sfd). 
Not only are the mean velocity profile and dif- 
fusion length reproduced, but so are various mi- 
croscopic statistics of the packing geometry, such 
as the two-body and three- body correlation func- 
tions l|Rvcroft et all l2005|) . These surprising re- 
sults are quite insensitive to the details of the 
relaxation step, apparently due to the very small 
and infrequent particle overlaps which arise from 
the cooperative mechanism in the spot simula- 
tion. Another, deeper reason may be that the 
geometry of dense flowing random packings has 
universal features, which are achieved by the to- 
tally different dynamics of spot and DEM simu- 
lations. 

4. Mathematical analysis of diffusion 

4.1. A stochastic differential equation 

In this section, we return to the general formu- 
lation of the Spot Model and analyze tracer dif- 
fusion in the continuum limit. It is clear from the 
simulations in Fig.|5]that the basic model in Fig-121 
gives a reasonable description of the dynamics of 
a single particle tracer, even though the multi- 
scale relaxation step in Fig. 0] is needed to pre- 
serve realistic packings. The relatively small size 
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Figure 5. Simulations of granular drainage in a quasi-two-dimensional silo. Top: T wo-dimensional 
simul ations using (a) the Spot Model without packing constraints and (b) the Void Model l|Bazant et all 
2005). Bottom: Three-dimensional simulations using (c) the Spot M odel with mult i scale relaxation and 
(d) the Discrete Element Method for frictional, visco-elastic spheres fo.vcroft et alll2005h . Particles arc 
colored according to their initial positions in horizontal stripes, lOd thick. 
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of the relaxation displacements makes it reason- 
able to regard them a small additional "noise" in 
a mathematical analysis of tracer diffusion. Here, 
we will neglect this small (but complicated) noise 
and view its average effect as incorporated statis- 
tically into the spot influence function, w(r p ,r a ), 
in Eq. igj. 

We begin by partitioning space as shown Fig.EI 
where the nth volume element, Avj n \ centered 
at Tg contains a random number, ANg n \ of 
spots at time t (typically one or zero). In a time 
interval, At, suppose that the jth spot in the 
nth volume element makes a random displace- 
ment, AR^'"' (which could be zero). According 
to Eq. the total displacement, AR p , of a par- 
ticle at v p in time At is then given by a sum of 
all the random displacements induced by nearby 
moving spots, 

AR p = -£ ]T W (r p ,rW+AR<^)ARM. 

Note that the spatio-temporal distribution of 
spots, ANg U \ is another source of randomness, 
in addition to the individual spot displacements, 
AR^' n ', so that each particle displacement is 
given by a random sum of random variables. 

In the limit of infinitesimal displacements, we 
arrive at a non-local, nonlinear stochastic differ- 
ential equation (SDE): 

dRp(t) = -JdN s (r s ,t) (15) 

w(r p (t),r s + dR s (r s ,t))dR s (r s ,t), 

where the stochastic integral is defined by the 
usual limit (infinitely refined partition of space) 
of the random Riemann sum in Eq. (|14|) . This 
equation diffe rs from standard nonlinear SDEs 
(Risken, 1996) in two basic ways: (i) The tracer 
trajectory, 

r p (t) = f dR p (T) (16) 

Jt=0 

is passively driven by a stochastic distribution 
of moving influences (spots), dN s (r s ,t), which 
evolves in time and space, rather than by some 




Figure 6. Sketch of a particle interacting with 
a collection of passing spots, showing some of 
the quantities involved in the nonlocal SDE, Eq. 

CSJ- 



internal source of independent noise, and (ii) the 
stochastic differential, dR p (t), is given by a non- 
local integral over other stochastic differentials, 
dH s (r s ,t), associated with these moving influ- 
ences, which lie at positions, r s , at finite distances 
away from the particle at r p . 

4.2. A Fokker- Planck equation 

In general, the various stochastic differentials 
in Equation l|15|l are correlated, which signifi- 
cantly complicates analysis. Here, we will make 
the reasonable first approximation of an ideal gas 
of spots, where the tracer particle sees an inde- 
pendent, random configuration of non-interacting 

spots at e ach infinitesimal time st ep. As in an 

ideal gas ([Hadiiconstantinou et all feOOffl , spots 
are thus distributed according to a Poisson pro- 
cess with a given mean density, p s (r s ,t). In ad- 
dition to neglecting correlations caused by inter- 
actions between spots, we disregard the follow- 
ing facts: (i) The distribution of spots in space, 
{dN s (r s , i)}, at time t depends explicitly on the 
distribution and displacements at the previous 
time, t — dt, via the spot random walks; (ii) Each 
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spot, due to its finite range of influence, affects 
the same particle for a finite period of time, so 
any persistence (autocorrelation) in the spot tra- 
jectory is transferred to the particles, in a nonlin- 
ear fashion controlled by w(r p , r s ). 

In the spot-gas approximation, the tracer par- 
ticle performs a random walk with independent 
(but non- identically distributed) displacements, 
which depend non-locally on a Poisson process 
for finding spots. Therefore, the propagator, 
P p (r, t\ro, to), which gives the probability density 
of finding the particle at r at time t after being at 
ro at time tn, satisfies th e following Fokker-Planck 
equation ijRiskenLlTflQ^ . 



BP, 

-flf + V-(u p P p ) 
with drift velocity, 

(dR P (r,t)) 



VV : (D P P P ), 



Up(r,t) 



dt 



lim 

At^O 



(ARp(r,t)) 
At 



(17) 



(18) 



and diffusivity tensor, 



Df (r,i) = 



2 eft 



(19) 

, ; — .) The 

Fokker-Planck coefficients may be calculated 
by taking the appropriate expectations using 



(Here VV : A denotes J2 a J2f 



d 2 A° 



in) 



and At 







Eq. in the limits AV S 
(in that order), which is straightforward since we 
assume that spots do not interact. Here, the spot 
displacements, AR^' (ri ), and the local num- 
bers of spots, ANs , are independent random 
variables in each time interval, and they are in- 
dependent of the same variables at earlier times. 

In order to calculate the drift velocity, we need 
only the mean spot density, p s (r s ,t), defined by 

The result, 



(AN^) = p s (r s ,t)AV s (n) . 



u p (r p ,t) = - J dV s w{v p ,v s ) [p s (r s ,t)u s (r s ,t) 
-2D a (r„f).Vp.(r a ,t)] (20) 

exhibits two sources of drift. The first term in 
the integrand is a particle drift velocity, which 
opposes the spot drift velocity, 



u s (r,t) 



(dR a (r,t)) 
dt 



(21) 



as in Eq. J3J. The second term, which depends 
on the spot diffusion tensor, 



Dl y, (r,t) = 



(dRV'dR 
2dT 



(22) 



is a " noise- induced drift" , typical of nonlinear 
SDEs ilRiskerl Hool . which causes particles to 
climb gradients in the spot density. This extra 
drift is crucial to ensure that particles eventually 
move toward the source of spots, e.g. the orifice 
in granular drainage. Both contributions to the 
drift velocity in Eq. (|20|l are averaged non-locally 
over a finite region, weighted by the spot influence 
function, w(r p ,r s ). 

In order to calculate the diffusivity tensor, we 
also need information about fluctuations in the 
spot density. From the spot-gas approximation, 
we have 

(AN^AN^) = S m , n ((AN^) 2 ) = 0((AV s ^) u ), 

where v = 1 for a Poisso n process and v < 1 for a 
hype r-uniform process l|Torauato and Stillingerl 
2003). It turns out that such fluctuations do not 
contribute to the diffusion tensor (in more than 
one dimension), and the result is 

T> p (r p ,t) = J dV s w(r p ,r s ) 2 p s (r s ,t)D s (r s ,t). 

(23) 

Note that the influence function, w, appears 
squared in Eq. (|23|) and linearly in Eq. I|20|) , which 
causes the Peclet number for tracer particles to 
be of order w smaller than that of spots (or free 
volume), as in Eq. (fTT|l . 

Higher-order terms a Kramers-Moyall expan- 
sion generalizing Eq. 1)17(1 for finite independent 
displacements, which do depend on fluctuations 
in the spot density, are straightforward to calcu- 
late, but beyond the scope of this paper. Such 
terms are usually ignored because, in spite of im- 
proving the approximation, they tend to produce 
small ne gative p r obab ilities in the tails of distri- 
butions (Risken, 199(J. In granular materials, 
however, velocity gradients can be highly local- 
ized, so the correction terms could be useful. 

4.3. Spatial velocity correlation tensor 

For any stochastic process representing the mo- 
tion of a single particle, it is well-known that 
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transport coefficients can be expressed in terms 
of temporal corr elation functi ons via the Green- 
Kubo relations For example, the 

diffusivity tensor in a uniform flow is given by 
the time integral of the velocity auto-corrclation 
tensor, 



dt{uz(t)u$(o)) 



(24) 



where U p (t) = {Up} — dRp/dt is the stochastic 
velocity of a particle. (A similar relation holds 
for spots.) 

In the Spot Model, nearby particles move coop- 
eratively, so the transport properties of the collec- 
tive system also depend on the two-point spatial 
velocity correlation tensor, 

_(U?(ri)C#(r a )> 

p 



a a/3 (ri,r 2 



(25) 



which is normalized so that C™^(r, r) = 1. We 
emphasize that the expectation above is condi- 
tional on finding two particles at ri and r 2 at a 
given moment in time and includes averaging over 
all possible spot distributions and displacements. 
Substituting the SDE JT2J into Eq. yields 

Cf(v u v 2 ) = f dV s p s (r s )w(r u r s )w(v 2 ,r s ) 

Dfiv^/^Df^Df^) (26) 

assuming independent spot displacements. 

Equation Q26|l is an integral relation for cooper- 
ative diffusion, which relates the spatial velocity 
correlation tensor to the spot (or free volume) dif- 
fusivity tensor via integrals of the spot influence 
function, w(r p ,r s ). If the statistical dynamics of 
spots is homogeneous (in particular, if D s is con- 
stant), then the relation simplifies: 

C^(ri,r 2 )= (27) 

J dv s p 3 (r s )w)(ri,r 3 ) 2 J dv^ p B {r' s )w(r 2 ,r' s ) 2 

if also the tensor is diagonal, Df" oc 6~ a ,p. If the 
statistical dynamics of particles is also homoge- 
neous, as in a uniform flow (p s = constant), then 
it simplifies even further: 

/ dV s w(r - r s )w(-r s ) 



Cf{r) = 



where we have assumed that the spot influence 
function, and thus the correlation tensor, is trans- 
lationally invariant (r = ri — r 2 ). In this limit, 
as mentioned above, the velocity correlation func- 
tion is simply given by the (normalized) overlap 
integral for spot influences separated by r. 

4.4. Relative diffusion of two tracers 

The spatial velocity correlation function affects 
many-body transport properties. For example, 
the relative displacement of two tracer particles, 
r = ri — r 2 , has an associated diffusivity tensor 
given by, 

D a P(r u r 2 ) = Df{v{) + Df (r 2 ) (29) 



2Cf{v ll v 2 )J Dr (v 1 )D^{v 2 ) 



In a uniform flow, the diagonal components take 
the simple form 

D aa (r) = 2 D^ a (1 - C% a {r)) (30) 

which may be used above to estimate the cage- 
breaking time, as the expected time for two par- 
ticles diffuse apart by more than one particle di- 
ameter. A more detailed calculation of the rela- 
tive propagator, P(r, t\r , to), neglecting tempo- 
ral correlations (as above) would start from the 
associated Fokker-Planck equation, 

dP(r, t) 



dt 



= VV: (D(r)P(r,*)) 



(31) 



/ dV s w(i 



(28) 



with a delta-function initial condition. (In a non- 
uniform flow, one must also account for noise- 
induced drift and motion of the the center of 
mass.) This analysis does not enforce packing 
constraints, so it allows for two particles to be 
separated by less than one diameter. A hard- 
sphere repulsion may be approximated by a re- 
flecting boundary condition at |r| = d when 
solving equations such as (|31|) . but there does 
not seem to be any simple way to enforce inter- 
particle forces exactly in the analysis. 

5. Tracer diffusion in granular drainage 

5.1. Statistical dynamics of spots 

The analysis in the previous section makes no 
assumptions about spots, other than the exis- 
tence of well-defined local mean density, mean ve- 
locity, and diffusion tensor, which may depend on 
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time and space. As such, the results may have rel- 
evance for a variety of dense disordered systems 
exhibiting cooperative diffusion (see below). In 
this section, we apply the model to the specific 
case of granular drainage, in which spots diffuse 
upward from a silo orifice, as in Fig. [5] Our goal 
here is simply to show how to derive continuum 
equations from the Spot Model in a particular 
case, but not to study any solutions in detail. 

For simplicity, let us assume that each spot un- 
dergoes mathematical Brownian motion with a 
vertical drift velocity, u s — v s z, and a diagonal 
diffusion tensor, 

























d! 



(32) 



which allows for a different diffusivity in 
the horizontal (_L) and vertical (||) direc- 
tions due to symmetry breaking by gravity. 
In that case, the propagator for a single 
"spot tracer", P s (x, z, i|xo, Zq, t), satisfies another 
Fokker-Planck equation, 



dt dz 



The coefficients may depend on space (e.g. larger 
velocity above the orifice than near the stagnant 
region), as suggested by the shape of s ome e xper- 
imental density waves l|Baxter et allflQSf^ . 

The geometrical spot propagator, 
"P s (x, \z, Xq, Zo), is the conditional probability 
of finding a spot at horizontal position x once it 
has risen to a height z from an initial position 
(x , Zq). For constant v s and D s , the geometrical 
propagator satisfies the diffusion equation, 



oz 



(34) 



where b — Dj:v s is the kinematic parameter. If 
spots move independently, this equation is also 
satisfied by the steady-state mean spot density, 
p s (x, z) , analogous to Eq. of the Kinematic 
Model. However, the mean particle velocity in 
the Spot Model, Eq. (|2*U)l . is somewhat different, 
as it involves nonlocal effects (see below). 



The time-dependent mean density of spots, 
p s (x, z,t), depends on the mean spot injection 
rate, Q(xo,ZQ,t) (number/areaxtime), which 
may vary in time and space due to complicated 
effects such as arching and jamming near the ori- 
fice. It is natural to assume that spots are in- 
jected at random points along the orifice (where 
they fit) according to a space-time Poisson pro- 
cess with mean rate, Q. In that case, if spots 
do not interact, the spatial distribution of spots 
within the silo at time t is also a Poisson process 
with mean density, 



din 



(35) 



p s (x,z,t) = / cZx / dz / < 
J J Jt <t 

Q(x , zo,<)P s (x, z,<|x , z ,t ). 

For a point-source of spots (i.e. an orifice roughly 
one spot wide) at the origin with flow rate, Qo(t) 
(number/time), this reduces to 



dt Q (to)P fl (x,z,t|0,0,to), 



t <t 



(36) 



where P s is the usual Gaussian propagator for 
Eq. (|33(l in the case of constant u s and D s . In 
(33) reality, spots should weakly interact, but the suc- 
cess of the Kinematic Model suggests that spots 
diffuse independen tly as a first appr oximation in 
granular drainage ( Choi et allEf 



5.2. Statistical dynamics of particles 

Integral formulae for the drift velocity and dif- 
fusivity tensor of a tracer particle may be ob- 
tained by substituting the spot density which 
solves Eq. ()35(l into the general expressions l|20|l 
and H2.'i[l . respectively. For example, if spots only 
diffuse horizontally (Z)J = 0), then the mean 
downward velocity of particles is given by 



v p (r,t) = J dV s w(r p ,r s ) p s (r s ,t)v s (r s ,t) 



(37) 



Note that the mean particle velocity is a nonlocal 
average of nearby spot drift velocities. 

For simplicity, let us consider a bulk region 
where the spot density varies on scales much 
larger than the spot size. In this limit, the in- 
tegrals over the spot influence function reduce to 
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the following "interaction volumes" : 
Vfc(r) = J dr s w{r,r s ) k 



(38) 



for k = 1,2. (Note that V 1 = V s above.) The 
equation for tracer-particle dynamics (|17|l then 
takes the form, 



dP p 

at 



d_ 

dz 



VsPs 



2D 



dp ( 
dz 



(39) 



-2V X • (D^(V ±Ps )V 1 P p ) 

(d\ Ps V 2 P p ) + Vl (D^p s V 2 P p ) . 

Again, it is clear that rescaling the spot density 
is equivalent to rescaling time. 

When the spot dynamics is homogeneous (i.e. 
u s and D s are constants), Equation Q39J1 simpli- 
fies further: 



1 dP p 

»v„ at 



d_ 

dz 



b^ 2 



■ bl 



i:) 2 



p dz 2 



(PsP P ) (40) 



-2b ± V-(P p V Ps )-2b^(p p ^ 
oz \ oz 

where b 1 - = b = Djr/v s and foil = d\/v s are the 
spot diffusion lengths and b p — b p V 2 /Vi and bp = 
b^Vz/Vi are the particle diffusion lengths. In this 
approximation, the latter are given by the simple 
formula, 



p 

b 1 - 



b p _ J dV s w(r,r s ) 2 



J dV s w(r,] 



(41) 



which generalizes Eq. Qllj l for a uniform spot with 
a sharp cutoff. The physical meaning of the dif- 
fusion lengths becomes more clear in the limit 
of uniform flow, p s = constant. In terms of the 
position in a frame moving with the mean flow, 



C = V 



where v p = v s V s p s , we arrive at a 



simple diffusion equation, 



(42) 



where £, the mean distance dropped, acts like 
ti me, consistent wit h the experimental findings 
of lChoi et all ll2004l) . 



6. Possible application to glasses 

We have seen that the Spot Model, in its sim- 
plest form, accurately reproduces the kinematics 
of bulk granular drainage, so it is tempting to 
speculate that it might be extended to flows in 
other amorphous materials. In this section, we 
briefly consider evidence for spot-like dynamics 
in glasses, but we leave further extensions of the 
Spot Model for future work. 

Experiments have revealed ample signs of "dy- 
namical heterogeneity" in supercooled liquids and 
glasses llHansen and McDonaldLlT^lKoTllT99l 
lAngell et all 1200(1) . but the direct observation of 
cooperative motion has been achieved only re- 
cently. Rat her than com pact regions of relax- 
ation, [DonatOi^n *|l998j) have observed "string- 
like" relaxation in molecular dynamics simula- 
tions of a Lennard-Jones model glass. The 
strength and length scale of correlations increases 
with decreasing temperature, consistent with the 
Adam-Gibbs hypothesis. Such cooperative mo- 
tion would be difficult to observe experimen- 
tally in a molecular glass, but lWeeks et all l|2000j) 
have used confocal microscopy to reveal three- 
dimensional clusters of faster-moving particles in 
a dense colloids. In the supercooled liquid phase, 
clusters of cooperative relaxation have widely 
varying sizes, which grow as the glass transition 
is approached. In the glass phase, the clusters are 
much smaller, on the order of ten particles, and 
do not produce significant rearrangements on ex- 
perimental time scales. 

These observations suggest that the Spot 
Model may have relevance for structural rear- 
rangements in simple glasses. String-like relax- 
ation is reminiscent of the trail of a spot, in Fig.|3J 
An atomically thin chain might result from the 
random walk of a spot, roughly one particle in 
size, but carrying less than one particle of free vol- 
ume. Larger regions of correlated motion might 
involve larger spots and/or collections of interact- 
ing spo ts. Some key featur es of the experimental 
data of IWeeks et all l|2000l) seem to support this 
idea: (i) Correlations take the form of "neigh- 
boring particles moving in parallel directions" , as 
in Fig. |2 and (ii) the large clusters of correlated 
motion tend to be fractals of dimension two, as 
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would be expected for the random- walk trail of a 
spot, as in Fig.Efb). For a complete theory of the 
glass transition, however, one would presumably 
have to consider interactions between spots and 
thermal activation of their creation, motion, and 
annihilation. 

7. Conclusion 

In this paper, we have introduced a mecha- 
nism for structural rearrangements of dense ran- 
dom packings, due to diffusing spots of free vol- 
ume. Even without inter-particle forces, the Spot 
Model gives a reasonable description of tracer dy- 
namics, which is trivial to simulate and amenable 
to mathematical analysis, starting from a nonlo- 
cal stochastic differential equation. With a simple 
relaxation step to enforce packing constraints, the 
Spot Model can efficiently produce very realistic 
flowing packings, as demonstrated by the case of 
granular drainage from a silo. The spot mecha- 
nism may also have relevance for glassy relaxation 
and other phenomena in amorphous materials. 

Regardless of various material-specific appli- 
cations, the ability to easily produce three- 
dimensional dense random packings is interest- 
ing in and of itself. Current state-of-the-art 
algorithms to generate dense random packings 
are artificial and comput ationally expens i ve, es- 
pecially near jamming llTora uato et all [2000: 
iKansal et all 120021: lO'Hern et all l2002l 12003^ 
A popular e xample is the molecular dynamic s 
algorithm of iLubachevskv and Stillingerl lfl990). 
which simulates a dilute system of interacting 
particles, whose size grows linearly in time until 
jamming occurs. For each random packing gen- 
erated, however, a separate molecular dynamics 
simulation must be performed. In contrast, the 
Spot Model produces a multitude of dense ran- 
dom packings (albeit with some correlations be- 
tween samples) from a single simulation, which is 
more efficient than molecular dynamics, since it 
does not require the mechanical relaxation of all 
particles at once. It would be interesting to char- 
acterize the types of dense packings generated by 
the Spot Model and compare with the results of 
other algorithms. 

The Spot Model (with relaxation) also provides 



a convenient paradigm for multiscalc modeling 
and simulation of amorphous materials, analo- 
gous to defect-based modeling of crystals. The 
iteration between global "mesoscopic" simulation 
of spots and local "microscopic" simulation of 
particles leads to a tremendous savings in com- 
putational effort, as long as the spot dynamics is 
physically realistic for a given system. For exam- 
ple, a simple extension of t he multiscale simula- 
tions of granular spheres bv lRvcroft et al. 1 11200.51) 
would be to different particle shapes, such as el- 
lipsoids, which have b een shown to pack m ore ef- 
ficiently than spheres IjDonev e t al., 2004a|). The 
only change in the simulation would be to modify 
the inter-particle forces in the relaxation step for 
a soft-core repulsion with a different shape. 

A more challenging and fruitful extension 
would be to incorporate mechanics into the mul- 
tiscalc simulation, beyond geometrical packing 
constraints. One way to do this may be use 
the information about inter-particle forces in the 
spot relaxation step to estimate local stresses, 
which could then affect the dynamics of spots. It 
may also be necessary to move particles directly 
in response to mechanical forces, in addition to 
the random cooperative displacements caused by 
spots. Such extensions seem necessary to describe 
forced shear flows in granular and glassy materi- 
als. For now, at least we have a reasonable model 
for the kinematics of random packings. 
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